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Abstract. Recently, an Enskog-type kinetic theory for Vicsek-type mod¬ 
els for self-propelled particles has been proposed [T. Ihle, Phys. Rev. 
E 83 , 030901 (2011)]. This theory is based on an exact eqnation for a 
Markov chain in phase space and is not limited to small density. Previ¬ 
ously, the hydrodynamic equations were derived from this theory and 
its transport coefficients were given in terms of infinite series. Here, 
I show that the transport coefficients take a simple form in the large 
density limit. This allows me to analytically evaluate the well-known 
density instability of the polarly ordered phase near the flocking thresh¬ 
old at moderate and large densities. The growth rate of a longitudinal 
perturbation is calculated and several scaling regimes, including three 
different power laws, are identified. It is shown that at large densities, 
the restabilization of the ordered phase at smaller noise is analytically 
accessible within the range of validity of the hydrodynamic theory. An¬ 
alytical predictions for the width of the unstable band, the maximum 
growth rate and for the wave number below which the instability occurs 
are given. In particular, the system size below which spatial perturba¬ 
tions of the homogeneous ordered state are stable is predicted to scale 
with %/m where M is the average number of collision partners. The 
typical time scale until the instability becomes visible is calculated and 
is proportional to M. 


1 Introduction 

The emergence of collective motion of living things, such as insects, birds, slime 
molds, bacteria and fish is a fascinating far-from-equilibrium phenomenon which has 
attracted a great deal of cross-disciplinary attention [T]. The study of these systems 
falls under the broader umbrella of “active matter”. This term stands for collections of 
agents that are able to extract and dissipate energy from their surrounding to produce 
systematic motion m- Non-living examples of active matter are chemically powered 
nanorods [1], networks of actin fibers driven by molecular motors and swarms 
of interacting robots m- In 1995, a minimal computational model that captures 
the essentials of collective motion without entering too much detail was introduced 
by Vicsek [5]. In this so-called Vicsek model (VM), agents are represented as point 
particles traveling at constant speed. The agents follow noisy interaction rules which 
aim to align the velocity vector of any given particle with its neighbors. Later, more 
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sophisticated models that include additional interactions such as attraction and short- 
range repulsion were introduced jhll 012111281 . On the theoretical side, coarse-grained 
descriptions of the dynamics in terms of phenomenological hydrodynamic equations 
have been proposed on the basis of symmetry and conservation law arguments.Well- 
knowns examples are the Toner-Tu theory for polar active matter mm and the 
theory by Kruse et al. for active polar gels m- While being very successful, a dis¬ 
advantage of these approaches is that no link between microscopic collision rules and 
the coefficients of the macroscopic equations is provided. In particular, since all coef¬ 
ficients must be related to only a few microscopic parameters, they cannot be varied 
independently and the actual parameter space should be more restricted than the 
hydrodynamic equations might suggest. Furthermore, by just postulating equations 
there is the immanent danger that some terms or even entire equations might have 
been ommitted, which could become relevant in new, previously untested, situations. 
To address the missing-link issue, several groups have derived hydrodynamic equa¬ 
tions from the underlying microscopic rules and provided expressions for the transport 
coefficients in terms of microscopic parameters |16I26I34I17I35I36I42] . One of the first 
attempts to put the Toner-Tu theory on a microscopic basis was published by Berlin 
et al., who studied a Vicsek-type model with continuous time-dynamics and binary 
interactions by means of a Boltzmann approach mm- While the authors recently 
clarified |23) that even at low density their underlying microscopic model is not iden¬ 
tical to the Vicsek-model, many predictions agree qualitatively with the ones for the 
VM. This Boltzmann approach was later extended to systems of self-propelled par¬ 
ticles with nematic and metric-free interactions, and hydrodynamic equations were 
drived |42I33] . Very recently, the very basis of these derivations - the Boltzmann 
equation and its validity in active matter - has been critically assessed [11120122] . In 
particular, it has been shown that, at least near the threshold to collective motion, 
the binary collision assumption is not valid in the VM at realistic particle speeds and 
even at very low densities [20] . Furthermore, it has been demonstrated that the mean- 
field assumption of molecular chaos is not justified near the threshold to collective 
motion in soft active colloids and that the Boltzmann theory must be amended by 
pre-collisional correlations m- A first attempt to rigorously include correlations by 
means of a ring-kinetic theory for Vicsek-type models was put forward by Chou et 
al. [31] . This theory is very complicated; work is still in progress to simplify it m- 

Nevertheless, arguments from ring-kinetic theory can be used to confirm the plau¬ 
sible presumption that mean-field theories become reliable in the VM at large parti¬ 
cle speeds (or time steps) and/or at large particle densities jU]. Here, large density 
means that the average number M of collision partners of a given agent is much 
larger than one. It is known from equilibrium statistical mechanics that the behavior 
of spin-models become mean-field or more mean-field-like if the number of interaction 
partners is increased which can either be achieved by extending the range of interac¬ 
tion or by increasing the spatial dimension. It appears reasonable to assume a similar 
tendency in the VM which consists of moving “spins”. 

The large density limit M ^ 1 is beyond the capability of Boltzmann approaches 
because Boltzmann equations are restricted to binary interactions. However, a re¬ 
cently proposed Enskog-like theory [171191201 has no limitation on density and has 
already be applied to a model with M = 2 ... 7 and metric-free interactions [53] . This 
Enskog-like theory is based on an exact equation for a Markov chain in phase space. 
In a previous paper hydrodynamic equations were derived from this theory and 
its transport coefficients were given in terms of infinite series. In this paper, I ana¬ 
lyze the transport coefficients in the large density limit, M ^ 1, and show that they 
take a simple form. This allows me to analytically evaluate the well-known density 
instability of the polarly ordered phase near the flocking threshold [30126] , and to ob¬ 
tain simple formulas and scaling laws for the dispersion relation. Note that a similar 
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analysis for the opposite limit M -C 1 has already been performed by Bertin et al. 


m- 


The large density expansion performed in this paper is supposed to be benefitial for 
several reasons. First, one does not have to worry about the validity of the underlying 
mean-field assumption of Molecular Chaos: As discussed above, this assumption is 
expected to be asymptotically valid for M —>• oo. Thus, the obtained relations should 
be quantitatively correct in this limit. As will be shown below, at large density, the 
restabilization of the homogeneous ordered phase happens closer to the threshold. 
Hence, this effect occurs within the validity domain of the hydrodynamic theory, and, 
aiming at quantitative agreement, one is not forced to evaluate the kinetic eqution 
directly as proposed in Ref. m- Second, biologically realistic models of fish schools 
m and bird flocks [35] assume that M is between 2 and 7, thus an approach for large 
M seems promising. Third, similar to the opposite limit, M —>■ 0, the large density 
approach leads to a simplification of the expressions for the transport coefficients and 
allows a physical interpretation of the instability close to the threshold of collective 
motion. 

An obvious concern about the large density limit is how realistic it is for active 
matter in general. Of course, for example, for systems of active colloids with short- 
ranged excluded volume interactions it will be problematic. However, it does make 
sense for systems with long-ranged alignment interactions. Speaking in terms of the 
terminology put forward by Couzin et al. |10j . the condition M 3> 1 can occur in a 
model with a large “zone of orientations” and a very small “zone of repulsion” Fur¬ 
thermore, in the process of doing the actual calculations and performing expansions 
in powers of 1/M, one realizes that error terms are often proportional to exp(—M). 
For example, for M = 2, this is already a reasonably small number. Therefore, it is 
possible that the large density formulas derived in this paper could still be valid at 
M slightly larger than one without loosing too much accuracy. 

To the best of my knowledge, almost no analytical results exist for the Vicsek- 
model at high particle density. This paper is intended to fill this void. One of the 
main results is the prediction that the homogeneous ordered phase remains stable to 
small perturbations if the linear system size is chosen to be smaller than L* = 'k\/2M. 
Another main result is that in systems larger than L* the number of time steps to 
develop an instability is equal to or larger than IGttM. 


2 Vicsek model 

In the two-dimensional Vicsek-model (VM) |8I40) . N particles with average number 
density po = N/A move at constant speed vq. The system of area A evolves in discrete 
time intervals of size r. The particles are assumed to have zero volume. The evolution 
takes place in of two steps: streaming and collision. During streaming, the particle 
positions Xj (t) are updated in parallel according to 


Xj (t + r) =Xj (t) -trvj (t). 


( 1 ) 


Because the particle speeds remain the same at all times, the velocities, Vj(t), are 
parametrized by the “flying” angles, 0j, Vj = Vo{cos6j,siii6j). In the collision step, 
the directions 9j are modified in the following way: a circle of radius R is drawn 
around the focal particle j, and the average direction (pj of motion of the particles 
within the circle is determined according to 



(2) 
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where the sum goes over all particles inside the corresponding circle (including particle 
j). The new flying angles are then given as, 

(^ + ''■) = (3) 

where is a random number uniformly distributed in the interval [— 77 / 2 , 77 / 2 ]. In this 
paper, a forward-upating rule [33], corresponding to the so-called standard Vicsek- 
model, is used: the updated positions Xj{t + t) instead of the “old” positions Xj{t) 
enter the calculations of the average directions <l>j. Important dimensionless parame¬ 
ters of the model are (i) the noise strength 77, (ii) the average number of particles in 
a collision circle, M = ttRI^pq, which is basically a normalized density, and (iii) the 
ratio 7 = R/X of the range R of the interaction to the mean displacement during 
streaming A = vqt which can be interpreted as a mean free path (mfp). 


3 Hydrodynamic theory 

3.1 Macroscopic equations for the Vicsek model 


Using Boltzmann’s assumption of molecular chaos, an Enskog-type kinetic equation 
for the one-particle density /(0,x,t) was derived from the exact evolution equation 
of the 7V-particle probability. This was done for the standard Vicsek-model (VM) [T7| 
as well as for other microscopic models [32124125] . By means of a Chapman-Enskog 
expansion and an angular Fourier transformation of the one-particle density, 

00 

f{e,x,t) = go{x,t) + ^ [gk{x,t) cos {k9) + hk{x,t) sin {k9)] (4) 

fe=i 

hydrodynamic equations were obtained for the dimensionless particle density p and 
momentum density w of the VM m- In this formulation, all times are scaled by the 
time step r, all lengths are scaled by the mean free path (mfp) , A = vqt, and all other 
quantities are scaled accordingly. For example, the actual two-dimensional number 
density p and momentum density w can be obtained from their scaled versions as 
p = p/X^ and w = wtiq/A^. To emphasize rotational invariance, the hydrodynamic 
equations are given in terms of the tensors H, Qi, and Q 2 which depend on spatial 
derivatives of p and w: 

dtp -I- V • w = 0 (5) 

dt'w -\-S/ ■ H = —bVp-\- (A — l)w -b Qi • w + Q 2 ■ Vp (6) 

with b = (3 — yl)/4. Here, A is the amplification factor for the first Fourier mode, 
which also determines the linear growth rate (1 — H) of the momentum density. This 
is because density and momentum density are the zeroth and first order moments 
of the one-particle distribution function, respectively, and therefore, the components 
of the momentum density are proportional to the first Fourier coefficients, Wx oc gi, 
Wy oc hi. In Ref. [HI, in the thermodynamic limit N ^ 00 , an infinite series was 
found for A, 
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dOn COS 6 cos 61 
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where 9 = 9{9i, 02,. ■ ■, On) is equal to the average angle defined in Eq. ([2]). Here, 
Mji is the local average number of particles within a collision circle centered around 
the position x, 


27r p 

ddf{9,x',t)= / dx'p{x.\t) (8) 

J |x—x^|<R 

For a homogeneous system, one has Mu = M = ttR^pq. If A is larger than one, the 
disordered state of zero average motion becomes unstable and a polarly ordered state 
with non-zero global momentum develops. Thus, setting A = 1 defines the threshold 
noise pciAI) of the flocking transition, at the level of homogeneous mean field theory. 
The hydrodynamic equations, Eqs. and are only valid in the vicinity of the 
transition to polar order, that is, for |A—1| <C 1. This condition was essential to obtain 
a closed set of just two equations, see Refs. [niD]. Physically, this assumption means 
that the first order Fourier coefficients gi and hi are evolving so slowly that the higher 
order coefficients, 92 , 53 , • ■ -, become enslaved to them [43) . Therefore, no additional 
equations for the temporal evolution of the higher order coefficients are needed. 

The momentum flux tensor H and the tensors Qi, <52, 




dx' 


Mx—x'|<R 


5 5 5 

H = ^hi fii Qi = ^ qi fli Qa = X! 

i—1 2 = 1 2 = 1 

are given in terms of five symmetric traceless tensors 

^a/3^7^7 

^2,oij3 — 

^3,oip — 

— '^a.^pP ’^p^aP P 

^5,aP = 2{dap){dpp) - Sapidjp)'^ . ( 10 ) 

The tensor f7i is the viscous stress tensor of a two-dimensional fluid. The transport 
coefficients in Eq. m were obtained in the additional limit of large mean free path, 
A = vqt ^ R, and are given in Table [TJ This means, contributions from the so-called 
collisional momentum transfer, which are supposed to be relevant at small mean free 
paths, see Ref. m, are neglected here. The rationale for this choice was that at small 
mean free paths and not too large density, the mean-field assumption (on which our 
theory is based) is supposed to break down anyway. 

The transport coefficients depend on the following variables, 


P 

q 

S 

r 


4 _ - 


-sin( 5 )^- ^Mn) 

7) n\ 

I — ^ 




-sin 

^ ' ':V 2 


\n-l)M^ '^J2{n) 


■ sin ■ 


E 


T] 2 ^ nl 

' n=2 


{n-l)M^-^Mn) 


871^7"^ •„ V ® 2 


„ ■ sin „ 

35 2 


E 

n—3 


(n- l)(n- 2 )M”-V 4 (n) 


( 11 ) 


where 7 is the ratio of the interaction radius to the mfp, 7 = R/tvq- Note that these 
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hj 


kj 

1 

i+p 

S 

s 

8(p-l) 

2(p-l) 

8(p-l) 

2 

p^-l-lOp-l-1 

S 

S(p-l-5) 

96(p-l)2 

4(p-l)2 

96(p—1)2 

3 

q 

p Sq 

P-1 

r Sq 

2(p-l) 

4 4(p-l) 

4 

4(p-l)2 

r Sq{p-3) 

2 2(p-l)2 

r Sq{p-4) 

12 12(p-1)2 

5 

q(p'^ + lOp+l) 
48(p-1)3 

r Sq(p^-2p+13) 

24 24(p-1)3 

Sq{p+5) 

48(p-1)3 


Table 1. The transport coefficients hj, qj and kj, defined in Eq. are expressed as 
functions of F, S, p, q, see Eq. CD- 

coefficients have a non-local dependence on position through the normalized density 
Mfi- The transport coefficients contain the following four types of n—dimensional 
integrals, 

1 /‘‘Z'K p‘2'K 

Jm{n) = /r, tn / d9i... dOn'I'm (12) 

Jo Jo 

where ifVn_ is given by S'! = cos^ 9 cos 29i, tp '2 = cos 9 sin 9 cos 9i sin 02 , 

if '3 = cos 9 cos 01 cos 202 , and >^'4 = cos 0 cos 0i cos 02 cos 03 . The average angle 0 = 

is a function of the angles 0i, 02 ,... 0„, see Eq. ©• 

The Navier-Stokes-type equation, Eq. (I 6 |) is consistent with the one from Toner- 
Tu theory mm but contains additional gradient terms. The additional terms are a 
result of a Chapman-Enskog expansion that includes all terms up to third order in 
a formal expansion parameter. Discussions of higher order expansions can be found 
in Refs. |2()l2ll22l2:f| . Eq. (|ni) has a simple homogeneous flocking solution: w = wo n 
and p = pq. The amplitude of the flow is given by 


Wo = 



(13) 


3.2 Large M expansion 


Using an analogy to the freely-jointed chain model for polymers, the angular integrals, 
Eqs. (na, can be asymptotially evaluated in the limit n ^ 1 [41]. To leading order, 
the results are m, 


Ji(n) 

J2{n) 


1 

8 n 

1 

8 n 


Jo{n) - 
Ji{n) ^ 


32n3/2 

30? 

32n3/2 


The angular integral in the definition of A, Eq. ([7]), is evaluated similarly. 


(14) 



I in) 


(15) 
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Even with these approximations, the infinite series in Eqs. © and (El), still cannot 
be calculated exactly. However, we note that they can be written as an average over a 
Poisson distribution. This is a consequence of the mean-field assumption of molecular 
chaos which predicts that the density fluctuations are equal to the ones of an ideal 
gas, e.g. are Poisson distributed. For example, using Eq. El we can rewrite the 
amplification coefficient A from Eq. © as 


/l = ^™(f)<(n+ !)'/»> 

where the average of an arbitrary function g(n) is defined as 


n—0 

OO 


n! 

rn— 1 


= E 


n—1 


-Mr 1) 

(n-1)! 


(16) 


(17) 

(18) 


Expanding g around the mean value M of the distribution one obtains, 

OO 

(5(n))=g(M) + ^gW(M)^ (19) 

k^l 

where g^^^M) is the k — th derivative of g evaluated at n = M, and /ifc are the 
so-called central moments, 

gk = {{n - Mf) (20) 

which are known for the Poisson distribution. For example, one has 


Ml = 0 

M2 = M3 = M 

fi4 = M + 3M^ ( 21 ) 


Thus, the following approximation is obtained, 

M 


{gin)) = g{M) -b —g''{M) 


-g"'{M) + ^ + 
6 ^ ^ ’ 4! 


■,(4) 


( 22 ) 


To evaluate transport coefficients, averages over polynomials in n are needed. There¬ 
fore, we investigate a generic example, 50 = {n + e)°‘ with some power a and a positive 
constant e. The approximation, Eq. (l2^ . then gives for small e M 


/ \^.n,ra(, , a(a-l) , a(Q;-l)(a-2) 0 ( 0 ; - l)(a - 2)(a - 3)[(1/M)-b 3] ^ 

{go)^M ^1 + —^ + - — -+--j 

(23) 

This approximation becomes accurate for M ^ 1. For very large M it suffices to 
replace {g) by g{M). Thus, replacing n by M in Eq. (ITCl) . leads to the leading order 
approximation for the amplification factor H, 


A 



V 



. 


(24) 
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Thus, fulfilling the critical noise condition T = 1 in the limit of infinite density 
M ^ oo requires that the critical noise rjc becomes equal to 27r. Expanding ijc 
around its infinite density limit gives, 


VC 



forM > 1 


Approximation (l25]l leads to 


(25) 



sin {r]c) 



(26) 


Introducing the relative noise distance i5 to the threshold. 


5 (27) 

VC 

we consider the difference of the amplification factor to its value at the threshold, 
A{r]c) = 1 and find. 


A-1 


271 . 


sin 


(l)-K¥l)«‘" 



for (5 <C 1. 


(28) 


As a preliminary to calculate transport coefficients, we express the variables p, q, S 
and r from Eqs. (HU, in terms of averages over the Poisson distribution. 


P 

q 

S 

r 


4 1 

-sin(p)((n + l)Ji(n + 1)) Ri — sin (p) 
T] 2r] 


47r'7^ 

-sin (r;)((n + 2 )J 2 (n + 2)) 

V 

sin ^((n + 2) J3(n + 2)) ? 

V 2 

gj sin |((n + 3) J 4 (n + 3)) 


— sin (rj) 

2r) 

7r3/2q2 ^ 

- 7= sin - 

Aqy/M 2 

7]-5/2^4 ^ 

-^ Sin — 

Arjy/M 2 


(29) 


and evaluate the averages to leading order by means of Eqs. da. Finally, using Eq. 
(l26ll . at the threshold one obtains 


p{vc) R 

1 

\/ttM 


q{vc) « 

2 / 71- 
^ \ M 


S{vc) « 

" ~4M 

(30) 

r{vc) « 

7r27^ 

4M 



These expressions, Eqs. (1501) are used to evaluate the transport coefficients from Table 
[T] in the limit of large M and at the flocking threshold. Only the leading order con¬ 
tributions in the limit of large M are given. For example, the product Sq ~ M“^/2 
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3 

hj 

(Ij 

kj 

1 

1 

8 

S 7r72 

2 8 M 

S vr7^ 

8 32M 

2 

1+lOp 1 

96 96 

S 71-7^ 

4 16M 

5 5.^, 57772 

96 384M 

3 

9 ~ 7^ 

2 2 V M 


r , Sq ^ 77274 

4 4 ~ 16M 

4 

1 ~ —'fi rn 

4 ^ 4 V M 

r 1 3Sq 7 r 274 

2^2^ 8M 

r , Sq 77274 

12 ^ 3 48M 

5 

g(l-l-lOp) ^ 72 yir 

48 48 V M 

r 1 13 Sg 

24 24 96M 

5Sq ^ 577^7® 

48 768M2 


Table 2. The transport coefficients hj, qj and kj, dehned in Eq. evaluated at the 
threshold rj = rjc, in the limit M ^ 1. 


in coefficients such as < 74 , fca and so on is neglected because F ^ M~^ decays less 
strongly in the infinite density limit. The results are given in Table [21 

It is instructive to evaluate the momentum density of a homogeneous, ordered 
system, close to the threshold, at d <C 1, According to Eqs. (fOl) . (pS)) . and Table |21 
we find 

Wo = (31) 


Even if the distance to the flocking threshold S, is tiny, the amplitude wq eventually 
diverges in the infinite density limit. It also has a strong dependence on the mean 
free path ratio 7 = R/X and diverges for infinite mean free path A. This seemingly 
unphysical behavior is just a consequence of the particular choice of dimensionless 
variables where all lengths are measured in units of the mean free path A, and all 
times in units of the time step r. The amplitude wq is a measure of global polar order 
but it is not a convenient variable in the large density limit. Instead, I “renormalize” 
Wq be relating it to the polar order parameter 17, which is defined as the average speed 
of a particle (|v|) divided by the individual particle speed vq. This way, 17 will always 
be between zero and one, and a small 17 <g; 1 is expected near the order-disorder 
transition. One has 


17 = 


< |v| > 

Vq 


Wq 

P 


(32) 


where p is the dimensionless particle density. It is related to the regular number 
density of a two-dimensional system, phy p = pX^ Using the definition, M = irR^p = 
ttR^p/X^ we find p = M/{'Ky‘^) and an expression for 17 is obtained, 


n 


Inserting wq from Eq. (EH), the order parameter near threshold follows as 


(33) 


17 = (34) 

In contrast to wq, the order parameter does not depend on particle speed. This is the 
expected result because, (i) we work in the mean-field approximation, and (ii) the 
collisional contribution to the transport coefficients was neglected earlier. 


3.3 Stability analysis 

A long wavelength instability has been reported in Vicsek-type models |2bl27ll7ll9| . 
This instability is strongest for longitudinal sound-wave-like perturbations and exists 
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in the ordered phase at a range of noise values rjs < rj < rjc, next to the flock¬ 
ing threshold. In this chapter, I reanalyse this instability in the large M-limit for a 
longitudinal perturbation with wave vector k = kfi. Without loss of generality, the 
collective motion is assumed to go into the x-direction, n = x. As small perturbation 
of amplitude Sp is imposed on the homogeneous particle density po- A similar per¬ 
turbation Sw is added to the x-component of the homogeneous momentum density 
Wx,o = where wq is given in Eqs. (US and (EH), 

p = Po + Sp e*'="+‘"* 

Wx=Wo+ Sw 6*'=“+'^* 

Wy = 0 , (35) 

and uj is the complex growth rate of the perturbation. From the continuity equation, 
Eq. a simple relation between Sp and Sw follows, 

Sp = ——Sw (36) 

UJ 

The hydrodynamic equation for the momentum density, Eq. (jH, in linear order of 
the perturbations leads to 

{uj — k'^hi + 2ikwQh^) Sw + {—ik^h 2 + ikw^h'^ — k'^w^hj) Sp = (37) 

{A — 1+ ikwoQi + Sw^qs) Sw + 

{—ikb + WqA' — k^WQq2 + ikwf^q^ -\- w^q'-^ + ikwgk^) Sp 


Here, all transport coefficients are evaluated for a homogeneous system with constant 
density po- The density derivatives of the coefficients A, and q^ are also needed 
and, for example, are given as 


A’ = 


dA 

dp 


Using Eq. El, one finds 


dA 


2r] 


2 ddi 

(38) 

M \2) 

(39) 


This expression is then evaluated at the threshold with the help of Eq. (El, leading 
to 




TT'J 


(40) 


Similiarly, one obtains 


q'sivc) 


dqs 

dp 


v=vc 


7r^7® 

, and h'^{qc) 


0 


(41) 


Substituting Eq. El into El, a quadratic equation for w is obtained, 

uj'^ + auj = fd (42) 


with complex coefficients /? = /3_r -|- i/S/, a = an + ioti- The real and imaginary parts 
of a and /3 are found as, 

q;_r = 1 — H — iw^qz — hik^ = —2wQq3 — hik^ 
a I = kwo{2h'i - qi) 

Pr = k^^—b u'q(p '4 -|- ^3 — /13) -|- fc^/12] 

Pi = kwQ[-A' - wlq'^ + k^{q2 - hj)] 


(43) 
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where Eq. (USD has been used to eliminate A. Eq. (HU) has two solutions, uj±, where 
UJ+ is the one with the larger real part. Solving for the real part, one finds. 


D / ^ an , (A^ + 

ReM = -- + -- 


1 + 


-T==) 

■</ + i/^ J 


1/2 


(44) 


with the abbreviations 


zi = 




4 

anai 


Pr 


-Pi 


(45) 


Instead of directly analysing this complicated expression, a different approach is pur¬ 
sued: Splitting the growth rate into its real and imaginary part, w = ujn + Eq. 
(gH) can be rewritten as two equations for uir and wj, 

^R~ A- ancoR — aiuj = Pr 
2ujruji - I - uruji + aiUjR = Pi ( 46 ) 


The sign of the real part lor determines whether the system is linearly stable. We 
therefore eliminate w/ from Eq. (|iS| and obtain a quartic equation with real coeffi¬ 
cients for the real part of w, 

040 ;^ -I- + C 2 a;^ + CiuiR = cq (47) 


with coefficients 


C 4 = 4 

C 3 = 8aR 

C2 = aj + 5a^ - APr 

Cl = aR[a^ + opR — APr\ 

Co = /5/ + ajanPi + PRa\ (48) 

This equation can be solved exactly by radicals as given in Eq. (l4^ . However, in the 
high density limit M 1 and near threshold 17 ^ 1, the terms in the quartic equa¬ 
tion can differ by several orders of magnitude and the analytical expressions become 
difficult to evaluate. Furthermore, our aim is to gain a better physical understanding 
and to obtain simple expressions for relevant quantities like the maximum growth 
rate of a perturbation of the ordered state. Therefore, in the next section, instead of 
relying on the exact solution of the quartic equation, various scaling regimes of ojr 
are identified by different Ansatzes to solve Eq. (HTl) approximately. 


3.4 Analysis of the dispersion relation 

3.4.1 The limit of vanishing wave number 

To see whether there is an instability at all, the limit of zero wave number fc —>■ 0 is 
investigated by inserting the “hydrodynamic” Ansatz 

UIR = dik^ + d2k‘^ + 0{k^) 


(49) 
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into Eq. dUl) and expanding the coefficients from Eq. (HRl) in powers of fc, 

C3 = C30 + 032*^ + 0{k‘^) 

C2 = C20 + C22fc^ + 0{k'^) 

Cl = Cio + Ci2fc^ + Ci4fc'‘ + 0(fc®) 

Co = co2fc^ + + coefc® + 0(fc®) (50) 

(51) 

with new coefficients that do not depend on k , for example, 

C 02 = wl {{A' + wlq'^f + 2wlq3{A' + wlq3){2h3 - qi) + Aw^qH-b + wl[q4 + ks - hij])} 
co4 = wl {2{A' + wlq3)(h4 - 52 ) + (gi - 2h3)[-hi{A' + Wogg) + 2wlq3[q2 - ^ 4 )] 

+453(10093/12 - bhi + tOo/ii [54 + k 3 - /13])} 

cio = -Swlql 

C 12 = wlqs {-2wl(2h3 - qif - 12wlq3hi + 8[-6 + 100(54 + fcs - ^ 3 )]} (52) 


Collecting terms of order k^ and k^, we find 


C02 

( 53 ) 

Cio 

— [C04C4Q — C02C20 — C10C12C02] 

Cio 

( 54 ) 


To simplify the analysis we evaluate the dispersion in the large density limit M ^ 1 
and replace ioq, using Eq. (IM)) by wg = f 7 M/( 7 r 7 ^). In order to trace the origin of 
stabilizing and destabilizing effects, mathematical “markers” are put on the transport 
coefficients. This is done by taking the asymptotic expressions from Table [5] and 
multiplying them with dummy variables that become equal to one in the limit M —> 
00 . For example, we have 


93 = 


/13 = — 


93 /^ 

' AM 

hf 


A' = 


93 = 


b = 


2y/^ 

2M 

q'gP 

8 M2 

b 


( 55 ) 


where / = and the “marker”-variables are given by 53 , / 13 , A' and so on. Inserting 
these “marked” coefficients into Eqs. (P|) and (ESI), one finds that for vanishing order, 
17 —>■ 0 , the coefficient di is always positive but becomes negative for a strongly 
ordered state, 17 « 1. To obtain quantitative results and to investigate the sign 
change of di, the dominant positive and negative terms must be balanced. Balance 
is achieved by assuming that the order parameter 17 is of order 1/M or smaller. 
Mathematically, this is imposed by introducing the scaled order parameter D = f2M 
and assuming that 17 is at most of order one. Replacing 17 = 17/M in the equations 
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for the coefficients Cj and Cij, Eqs. (PI) and dsa), and neglecting terms of higher order 
in 1/M, the coefficients dramatically simplify to, 


C02 

C04 

C06 

ClO 

Cl2 

Ci4 

C20 

C22 

C30 

C32 


a^qsbhi 
~ __16M 
bh\ 

8M3_ 

_M 

bhi 

~T 

4M2 

2b 

417^53 

M 

hi 


(56) 


Note that the mean free path dependence through the factor / = Try^ dropped out 
exactly from these expressions at any density. By means of Eqs. (1531) and ((Ml) the 
low wave number approximation, Eq. (03), can now be easily evaluated with the 
coefficients di and d 2 found as, 


di 

d2 


MB 


M^B 


n^ql 


B = 2{A'f 



bqlQ^ 


with 


(57) 

(58) 


While 6 , 53 and A' are just marker variables that are equal to one, their appearance 
tells us that the low wave number behavior fc —> 0 is controlled by only four transport 
coefficients: b, 53 , yl' and A. The latter coefficient enters implicitly because according 
to Eq. (TT^ . it is involved in the calculation of the order parameter 17. The most 
interesting prediction of Eq. dSZl) is that there is a critical order parameter 


nc = = ^2 (59) 

Vbqs 

below which di is positive. Thus, very close to threshold, for 0 < 17 < v^/M, we do 
have a long wave length instability but there is restabilization slightly further away 
from threshold, at 17 > '/2/M. Therefore, the width of the instability window (in 
12-space) scales as 1/M. Given that 17 is at most of order one, one sees that inside 
the instability region, the growth rate w increases rapidly with M since di ~ M and 
d 2 ^ m 3. In the high density limit M —>■ 00 this suggests that the point k = 0 
becomes singular, e.g. the value for uir at k slightly above zero is very different from 
uj{k = 0) = 0. This behavior is supported by the direct solution of Eq. (iTTl) . Moreover, 
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inside the instability window (where 17 < v^), the small wave number expansion is 
only valid for <C fci where ki is defined through the equality of the first two terms 
in the expansion of ojr, \di\k\ = 1 ^ 21 ^ 1 - This gives, 



Evaluating fci near the threshold to collective motion, that is for 17 <C 1, we have 
B ft! 2{A')‘^ and fci follows as 



(61) 


which goes to zero in the limit M —?> oo and also vanishes for 17 —>■ 0. Direct solution of 
the quartic equation (H7)) shows that the low k expansion of ujr is of very limited use 
because fci goes to zero too rapidly. In particular, the shape of the dispersion relation 
near the onset of collective motion, which is characterized by the maximum (positive) 
growth rate uin^max and the highest unstable wave number fco (where w/j(/co) = 0 ) 
cannot be determined using a Taylor-expansion of (jJr. This is because fco will be 
shown to scale with a smaller power of QjM. Therefore, in the large M limit, fco is 
outside the range of validity of the Taylor series, that is fco ^ fci. Different techniques 
to find /cq are presented in the next section. 

Furthermore, Eq. (1571) also tells us that the mathematical reason for the instability 
is a nonzero density derivative A' and that the sign of A' does not matter. One also 
sees from Eqs. dSZl) and (l 58 ll that it is the combination of the coefficients b and 
53 that initiate restabilization. The coefficient b describes the pressure term in the 
Navier-Stokes-like hydrodynamic equation ([5]) and is equal to the square of the speed 
of sound, whereas 53 controls the nonlinear (cubic) stabilization of the momentum 
density. 

3.4.2 The shape of the dispersion relation at large wavenumbers 

In the previous section we discussed that the Ansatz ujfi = + ^ 2 ^^ + • ■ • only 

gives a condition at what values of the order parameter 17 an instability occurs but 
fails to predict the shape of the dispersion relation inside the instability region. In 
this section, different scaling analyses of the quartic equation (I47p are presented that 
are valid at larger wave numbers k. In particular, k is assumed to be larger than 
the limit value ki from Eq. (EH), and hnally large enough to describe the value fcp 
above which ujr becomes negative again. This value is of particular interest because 
through the length scale Tmsi = 27r/fco it gives us an crude idea above which system 
size L* the high density Vicsek model developes soliton-like density waves [snnn]. 
The maximum value of the dispersion relation, ujR^max, is important to provide a 
lower bound ^TrliUR^rnax on the time scale such waves can pile up. 

To calculate fco and uiR^max the following scaling Ansatz is made, 



(62) 


(63) 


which will be justified a posteriori. The rescaled wave number k and the rescaled 
growth rate 6j are assumed to be of order one. Inserting these definitions into the 
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quartic equation dUl) and only keeping terms in the lowest order of the small param¬ 
eter 1/M, a simple quadratic equation emerges, 


UJ 


■ UJ 


172^3 




Q'^B uhl 

166 32 256 ' 


Its solution, written in terms of the unsealed variables, is 


(64) 




fi |yl'| 172^3 

mTH ~ 4M 


4M 


[17(7 — 


k^hi 

k^hi 

16 


(65) 


where we used the definition of the critical order parameter, Eq. (I59|) in the second 
line to visualize that lor can only be positive if 0 < 17 < I7c- This expression gives 
as a first estimate of the maximum value of uir, 


^R^max 


f? 

4M 


V2-f2 


( 66 ) 


and justifies tlie initial scaling choice, ujr ~ 1/M of Eq. (15^ . Very close to the 
threshold, at 17 <C 1, one finds. 


n _ f2 
VsM ~ 7i 


(67) 


with 12 = Mf2 and I7c = \/2. Clearly, expression (1651) ceases to be valid at very small 
k because in contrast to Eq. (091), it does not predict ujR{k = 0) = 0. Comparing 
neglected terms in Eq. (H7)) with the surviving ones in (1641) gives the condition 


fc » 6:2 


17 


— for 17 <C 1 


( 68 ) 


for approximation (l65l) to be valid, assuming proxirnity to the threshold to ordered 
motion. Near the restabilization point where 17 Ri Dq the wave number restriction 
changes to 


k 


/X 

M V 263 


(69) 


which can be realized by noting that 


B = hql{f2l - Q^) Ki 2f2chql{f2c - 12) (70) 


and assuming (l7c — 17) <C 1 when analyzing the terms in Eq. (ITTl) . 

Einally, setting w/j = 0 in Eq. (05]) leads to an estimate for k^, the wave number 
which delimitates the domain of unstable modes, 


fco 


/4Cg3[l7c - 17] 4121^2-12] 


Mhi 


= 2 \ 


M 


which for 17 -C 1 gives. 


ko 


V M /iy2^i/4 



(71) 


(72) 
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The result fco I/VM confirms the scaling Ansatz, Eq. (IMll . It is also in the range 
of validity, /cq 2> ^ 2 ) of the approximation (16511 because 17/M is always much smaller 
than a/ D/m for M ^ 1 and 17 < \/2. 

Maximizing fco and ujuik^) at fixed M but variable 17 leads to the “most danger¬ 
ous” value fio = •^c/2 = According to Eq. (IM)) this corresponds to a relative 
noise distance So, 


Sd = <5(f7£)) = 

from threshold. Eor example, at M = 5, the instability is strongest at 5 = i5ir 
Furthermore, calculating 


4\/4M^ 


(73) 

0.013. 


koi^o) = 

^maxi^o) = (74) 

we get a lower bound for the system size L* above which the homeogeneous ordered 
state is unstable and we can extract a typical time scale for the formation of this 
instability, T*, 

L* ^ tvq = tvottVqM (75) 

ko{s2j:,) 

Stt 

T* Ri T- = IGttM t (76) 

Here, we had to reinsert the time step r and the particle speed vq because fc and uj are 
dimensionless quantities. Thus we arrive at the prediction that the lower bound for 
L* increases proportional to v^M. Comparing the prediction for L* from Eq. (1751) to 
the insert in Fig. 1 of Ref. m very good agreement with the lower curve is found for 
M >2. This figure was obtained directly from the full dispersion relation, evaluated 
by Mathematica, and without making any large density approximations. We also find 
that the number of iterations to see the instability (in a system that is significantly 
larger than L*), is proportional to M. 


3.4.3 The dispersion relation at intermediate wavenumbers 


So far, we have found an approximation for loh given by Eq. (03) for small fc <C fci 
and another one, Eq. (163, for large k ^ k^. However, between fci and fc 2 there is a 
large gap in wave number space. The closer one is to the threshold value, the more 
this gap widens up since the ratio fc 2 /fci = diverges for 17 0. The previous 

approximations did not allow us to find the most unstable fcmax- One of the goals 
of this section is to find kmax which is expected to be inside the unexplored gap in 
wavenumber space. 

A direct evaluation of the quartic equation (H71) using Mathematica indicated that 
the cubic and the linear term do not seem to be relevant in this intermediate range 
of wave numbers, fci <g; fc <g; fc 2 . To put this observation on a more solid ground, the 
following scaling Ansatz is made 


UJ 

k 

Mi+“ 



(77) 
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with unknown positive exponents a, P and (j) and inserted into Eq. dSD. To obtain 
a meaningful balancing of terms in this equation and to prevent contradictions the 
coefficients must fulfill the conditions, 

1 Ot 

(/)=-(a + /3) and — < j3 < a (78) 

Z O 

At the lowest order in 1/M one hnds a balance between the quartic term ^ and 
the 0(a;°) term, leading to the approximation 


UJR K. 




1/4 




(79) 


Near threshold, 17 1 and setting the dummy variables to one, this becomes 


= (80) 

Away from the restabilization region, that is assuming B = 0(1) and 17 < 1, the 
approximation (1791) is valid in the following window: 


M 


/32 , 17 



(81) 


For 17 <C 1 and setting the dummy variables to one this window becomes 


, 4173 , 


(82) 


This means the scaling regime lor ^ \/k is only clearly visible very close to the 
threshold where 17 is smaller or equal to about 0.1. Note, that ki happens to be equal 
to fe 2 , and ks has the same scaling ~ /M than ki. 

The approximation for ujr at intermediate wavenumber, Eq. o , increases mono- 
tonically with A:, whereas the expression for high A:, Eq. dMl), decreases monotonically. 
However, the regions of validity for these approximations do not overlap. Therefore, 
they should not be equated to obtain an estimate of k^ax where ujr reaches its max¬ 
imum. To estimate kmax, we first construct another approximation whose region of 
validity extends above k^. This is done by not only balancing terms of order and 
in Eq. (ITTI) . but also including a contribution proportional to lo\. This leads to a 
quadratic equation in 


with the solution 

+ + (84) 

Compared to Eq. dSH), the range of validity is extended to 



fca < fc < 4^/3 


= ^5 forl7 <C 1 


(85) 






18 


Will be inserted by the editor 


Now, we have an overlap with the large k approximation since k 2 f2/M <C A :5 ^ 
(17To match the results in the overlap region, the square root in Eq. (151)) is 
expanded for large k with the result: 




8M^b 




0{l/k^) ioT^J^<^k^k5 


From the high k expression (1651) one finds 

2 _ \A'\^f2'^ k'^hi\A'\D 

8M% 8MV% 


for k ^ k 2 


( 86 ) 


(87) 


where a term proportional to k'^ could be neglected for k ^ k^ax ■ The two expressions 
dHU) and dlzl) becomes equal at the particular wavenumber kmatch, given by 


kmatch 


1 



( 88 ) 


It seems plausible to assume that kmatch is a good estimate of the wavenumber kmax- 
This means kmax is expected to be proportional to (17/M)^/^. This is confirmed by 
a more accurate calculation of kmax and ujmax for 17 <C 1: 


1 17 1 17 

“ ~ 31/227/4 M 


(89) 

(90) 


The prefactor in Eq. (1551) is about 10% smaller than the one in Eq. (1551) but the 
scaling is the same. 


3.4.4 Verification of scaling predictions 

To verify the proposed scaling laws, Eq. (H^ is evaluated numerically. Using Eq. (IT5)) 
the threshold value rjc for a given M is found from the condition A{rjc) = 1- For a 
given relative distance 5, the noise is determined by 77 = riciA — ^)- The transport 
coefficients are obtained from Table [TJ by inserting the large M approximations for 
p, q, r and S from Eqs. (ESI), into the corresponding expressions. Note that the 
approximations from Eq. (l30l) and from Tableware not used. Fig. [T] shows the scaled 
real part of the growth rate, = Wfl-\/8/17, as a function of scaled wavenumber 
ks = fc2“®/^l7“i/2 for M = 7 and i5 = 10“®. The values fco = 0.15566, kmax = 
0.01372 and ujR^max = 0.001233 can be read off the plot. Eqs. (1551) and (IT^ give 
the order parameters 17 = 0.004018 and 17 = 17M = 0.02813 which are used to 
calculate the limiting wavenumbers, ^2 to k^. Following Eqs. dMl), ESI) and dSSl) one 
finds k 2 = ki = 0 . 02 , k^ = 1.27 x 10“®, and k^ = 0.0482. As expected, the conditions 
ko ^ k 2 , and k^ kmax "C k^ are fulfilled. The limit fci = 7.1 x 10“7 is too small 
to play any role in the plot. For verification, Eqs. (ED, Ell) and EQl) are evaluated 
and give ko = 0.1507, kmax = 0.01322 and WR^max = 0.001377. These values show 
excellent agreement with the ones from Fig. El 

1 also checked a denser system with M = 100 and 5 = 10“® corresponding to 
17 = 0.05741. Again, the errors of the predictions for ko, kmax and uJmax are only 
between 2 and 5 percent. Finally, a very dense system with M = 8000, 6 = 10“® and 
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Fig. 1 . The scaled growth rate lur^s = u!r\/8/0 is plotted versus the scaled wavenumber 
ks = for M = 7 and M = 100. UJR was obtained from direct solution of Eq. 

The large k approximation, Eq. IMIl, is given by the dashed line and represents an 
inverted parabola. The scaling was chosen to test the predictions for ko and lor from Eqs. 
([72|l and (|67l). 


= 0.304 was tested. Here, the error of the prediction for was only 0.6%. Since the 
scaling exponent of 3/4 in the prediction for kmax seems unusual, I checked whether 
the exponents 1/2 and 1 could also make sense. It turns out that they can be ruled 
out because for the M = 7 sample, kmax would change by a factor of ~ 4=*=^, 
and for the M = 100 system kmax would change by a factor of ~ 6.5^^. These factors 
are in clear contradiction with the small errors of the kmax oc prediction. 


3.5 Discussion of the dispersion relation 


Close to the threshold to collective motion, that is for 17 = QM <C 1, the dispersion 
relation has a range of unstable modes. One of the main results of this paper can be 
seen in Fig. [TJ For M 1, the unstable modes are asymptotically described by an 
inverted parabola for the real part of the growth rate ojr. It is not a perfect parabola 
since close to zero wavenumber, wr abruptly turns down. Inside this region of rapid 
change, a power law applies, wr ~ The imperfection of the parabola becomes 
smaller and smaller with increasing density, as shown in Fig. [T] Loosely speaking, I 
found a “boundary layer” in the dispersion relation. For extremely small wavenumbers 
k K, Q there is another power law regime uj ^ k'^. However, for ilM ^ 1 and M ^ 1 
this regime is basically irrelevant. For example, even for M = 7 it is not visible in 
Fig.ffl 

Another main result is the scaling of the most unstable mode kmax ^ 17^/^ and 
^R,max 12, as well as the scaling of the range of unstable modes, feg ~ 17^/^. The 
different scaling laws for ko and kmax-, near the threshold to collective motion can be 
combined in the relation 



h.2 

^max 


= const. 


(91) 


Note, that this disagrees with the low density Boltzmann approach for a Vicsek-like 
model by Bertin et al [27]. The scaling reported there corresponds to k^/kmax = const. 
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4 Conclusion 

Previously, an Enskog-type kinetic theory for Vicsek-type models was proposed and 
hydrodynamic equations were drived HZ]. The corresponding transport coefficients 
were given in terms of infinite series which are difficult to handle. The fact that 
this theory is not restricted to low density has not been fully utilized yet. In this 
paper, I exploit the special properties of the Poisson distribution to obtain simple 
approximations of the transport coefficients. These expressions become exact in the 
infinite density limit but, in practice, are expected to be quite good as long as M, 
the average number of collision partners, is of order one or larger. Analyzing the 
hydrodynamic theory of the VM in the large M-limit is not only advantageous with 
respect to having simpler expressions but also because the underlying mean-field 
assumption is supposed to become exact near the flocking threshold for M —>• oo. 
Ranking transport coefficients by powers of the small parameter 1/M allows me to 
analytically evaluate the well-known density instability of the polarly ordered phase 
near the flocking threshold. 

The growth rate w of a longitudinal perturbation is calculated and a band of 
unstable modes with wavenumbers 0 < fc < fco is found. Some of the main results 
of this paper are given in Eqs. dzu), (ESI) dsni) which describe the scaling behavior 
of fco as well as the maximum growth rate and the most unstable mode number 
kmax in terms of density and size of the order parameter Q. It is found that there 
is only an instability for 0 < 17 < \/2/M. This means, for large M, a restabilization 
of the ordered phase occurs very close to the threshold - the size of the instability 
window in 17-space shrinks as 1/M. Thus, one can assume that for large enough 
M, the restabilization can be described within the validity domain of hydrodynamic 
theory. Inside the instability window, the largest value for fco is found for an order 
parameter value 17 = l/(M-\/2). This allows the calculation of the (approximate) 
maximum system size L* ~ below which the homogeneously ordered phase is 
stable. A corresponding time scale for the formation of density inhomogeneities was 
also calculated and found to increase linearly with density. The estimate for L* is in 
agreement with an earlier numerical evaluation of the dispersion relation without the 
high density approximation, see Ref. m- 

Furthermore, I show that the real part of the growth rate follows three different 
power laws: at very low wavenumber fc, one has i?e(w) ~ followed by a power 
law regime with ~ v^. At large enough wavenumber, one finds Re{w) ^ const — 
which represents an inverted parabola. The closer the system is moved towards the 
threshold, the more room in fc-space is taken by the inverted parabola, at the expense 
of the other two regimes. It is also interesting to see that the key features of the 
dispersion relation are determined by only five ( 6 , < 73 , /ii. A, and |yl'|) of the many 
transport coefficients of the hydrodynamic equations. 

Support from the National Science Foundation under grant No. DMR-0706017 is 
gratefully acknowledged. 
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